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ABSTRACT 


»  -  •] 


A  theoretical  model  capable  of  describing  the  flutter  and  divergence 
of  a  tube  conveying  high  velocity  fluid  is  developed  in  this  study.  The 
nonlinear  effect  due  to  tension  induced  by  bending  is  considered  along  with 
other  linear  forces.  The  partial  differential  equations  of  motion  are 
reduced  to  a  set  of  ordinary  differential  equations  by  Galerkin's  method  and 
the  time  histories  of  the  tube  motion  are  evaluated  numerically.  Numerical 
results  for  the  thresholds  of  flutter  instability  and  the  mode  shapes  of 
a  cantilevered  tube  are  in  good  agreement  with  those  obtained  theoretically 
and  experimentally  by  other  investigators.  (1) 

A  static  or  divergence  instability  occurs  at  a  certain  critical  fluid 
velocity  if  the  tube  is  simply  supported  at  both  ends.  At  higher  fluid 
velocity,  the  tube  becomes  dynamically  unstable.  The  limit  cycle 
oscillations  of  the  cantilever  tube  and  the  characteristics  of  the  simply 
supported  tube  motion  predicted  by  the  present  nonlinear  analysis  remain  to 
be  verified  quantitatively  by  experiment.  Plans  have  been  made  to  set  up 


an  appropriate  experiment. 
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NOMENCLATURE 


ajk  mode  shape  integral  given  by  Eq.  41 

a  cross  sectional  area  of  tube 

A  notation  for  matrix 

b^  mode  shape  integral  given  by  Eq.  42 

chl-  mode  shape  integral  given  by  Eq.  A-15 

ds  differential  length 

D  damping  factor  given  in  Eq.  9 

E  Young's  modulus 

fdamp  damping  force 

g  acceleration  due  to  gravity 

I  moment  of  inertia  of  tube  section 

i  length  of  tube 

m  mass  per  unit  length  of  tube 

mf  mass  per  unit  length  of  fluid 

q  column  vector  of  q^ 

gk 


generalized  coordinate 


some  constant  q^ 


radius  of  gyration 
end  reaction  force 
radius  of  tube 
inside  radius  of  tube 
outside  radius  of  tube 
time 

thickness  (Appendix  III) 
kinetic  energy 
tension 

potential  energy 
velocity  of  fluid  in  tube 
work 

coordinate 

displacement  of  tube 


nondimensional  velocity  = 
strain  (x-direction) 


3.! 


kth  cantilever  beam  mode  given  in  Eq.  34 


vi 


n^*1  simply  supported  beam  mode  given  in  Eq.  63 


variation 


delta  function 


operator,  see  Eq.  55 


critical  damping  ratio  of  itn  mode  without  fluid 


mode  shape  constant  of  Eq.  35 


mass  ratio  = 


m^  +  m 


viscous  damping  parameter  =  7.04  ^  /1-y 


nondimensional  coordinate  =  x/l 


nondimensional  displacement  *  y/r 


mode  shape  constant  of  Eq.  35 


stress  (x  direction)  in  psi 


nondimensional  time  =  / - — — -r  t 

1/  (m-+m)jr 


9-j,  &2  angles  of  deflection 


frequency 


natural  frequency  of  icn  mode  without  fluid 


nondimensional  frequency  = 


r(mf+mH 


3.52  u 


2 


As  early  as  1950,  the  observation  of  bending  vibration  of  a  simply 
supported  petroleum  pipeline  prompted  scientific  investigations  (2).  More 
recently,  it  was  observed*  that  a  short  pipeline  vibrates  in  its  second 
circumferential  mode  at  fairly  high  frequency  (300-800  Hz)  above  a  certain 
critical  flow  velocity  (see  Figure  1).  The  vibrations  may  sometimes  generate 
a  shrill  sound  (3).  If  the  pipeline  is  sufficiently  long,  this  instability 
is  then  superposed  on  the  flexural  oscillatory  instabilities  (Figure  2). 

High  performance  liquid-propellant  launch  vehicles  may  require  rapid  trans¬ 
fer  of  large  quantities  of  fluid  from  tanks  to  pumps  through  pipes  with 
relatively  thin  walls.  Thus,  a  study  of  the  influence  of  fluid  velocity 
on  the  static  and  dynamic  stabilities  of  propellant  lines  (tubes)  is  necessary 
(4).  Other  applications  occur  in  hydraulic  lines  and  human  lung  airways  (5). 

The  purpose  of  this  paper  is  to  analyze  the  observed  motion  of  an 
elastic  tube  conveying  fluid  using  a  mathematical  model  based  upon  the 
following  assumptions: 

(a)  The  three  basic  classes  of  operative  forces  are  inertial,  elastic 
and  aerodynamic  forces  coupled  by  the  elastic  deformations  of  the 
tube. 

(b)  The  fluid  being  conveyed  in  the  tube  is  inviscid,  incompressible 
and  non-heat  conducting. 

(c)  The  tube  may  be  approximated  by  a  uniform  beam. 

Earlier  investigations  have  developed  linearized  mathematical  models 
capable  of  describing  the  behavior  of  the  tube  up  to  and  including  the 
threshold  of  instability.  Section  1  of  Part  II  briefly  discusses  the  approach 
taken  by  Alan  S.  Greenwald  and  John  Ougundji  (1).  To  predict  the  large 
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amplitude  periodic  motion  after  the  threshold  of  instability  is  exceeded, 
nonlinearities  have  to  be  taken  into  consideration.  One  such  nonlinearity 
is  due  to  tension  induced  by  bending  and  consequent  stretching  of  the  line. 
Its  existence  is  closely  examined  in  Part  II,  Section  1.2.  So  far,  this 
study  is  believed  to  be  the  only  one  conducted  on  the  nonlinear  tension 
effect. 

Computer  simulation  is  used  to  solve  the  set  of  ordinary  differential 
equations  obtained  by  applying  Galerkin's  method  to  the  partial  differential 
equations  of  motion.  Data  have  been  analyzed  and  the  results  are  presented 
graphically  for  two  sets  of  tube  boundary  conditions,  cantilevered  (clamped- 
free)  and  simply  supported-simply  supported.  The  clamped-free  tube  under¬ 
goes  limit-cycle  oscillations  of  constant  peak  amplitude  and  frequency  at  a 
fixed  fluid  velocity.  Unlike  what  was  expected,  the  maximum  stress  does 
not  always  develop  at  the  root  of  the  cantilevered  tube  (the  clamped  end). 

In  fact,  the  point  of  maximum  stress  travels  along  the  length  of  the  tube 
during  one  period  of  the  limit  cycle  oscillation.  The  stress  and  deflection 
of  the  tube  increases  as  fluid  velocity  increases.  For  a  simply  supported- 
simply  supported  tube,  a  static  deflection  occurs  above  the  threshold  of 
instability  whose  amplitude  increases  as  the  flow  velocity  increases.  At 
still  higher  velocity,  the  tube  enters  into  a  dynamically  unstable  region. 


Section  1 


Derivation  of  Equations  of  Motion 


1 .1  General  Approach* 


To  obtain  a  differential  equation  of  motion  describing  the  tube,  an 
energy  approach  has  been  utilized  to  account  for  all  the  major  forces 
acting  on  the  nonconservative  system.  Hamilton's  principle  states 
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(T-U)  dt  + 


6W  dt  =  0 


(1} 


The  total  kinetic  energy  of  the  system  is  the  sum  of  the  kinetic  energy  due 
to  the  motion  of  the  line,  and  the  kinetic  energy  due  to  the  flowing  fluid. 
In  other  words. 


^  -  Tline  +  ^fluid 

(2) 

where 

£ 

Tline  =  7/  m(!t)2  dx’ 

(3) 

Tfluid  ‘7/  mf|7l2  dx 

14) 

In  this  paper,  the  x  and  y  components  of  the  fluid  velocity  are  approximated 
as** 


Vxn,V 


V  *&  +  V  &- 

y  at  9x 


(5) 


*  A  similar  discussion  was  previously  reported  in  reference  1. 

**  See  figure  3. 
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Then 


(6) 


o 


r- 


K.  • 


1  ^ 

Tfluid  =  7  /  m 

and  T  r  Wff)2  +  mf[V2  +  (§*)2  ♦  2V  §£f£  *  V2(|^)2]}  dx  (7) 

The  total  potential  energy  of  the  system  consists  of  elastic  strain 
energy  of  bending  and  the  elastic  energy  stored  in  tension  due  to  gravity. 

By  assumption  (c)  of  Part  I,  the  concept  of  simple  beam  theory  may  be  applied 
In  mathematical  notation  the  total  potential  energy  may  be  written  as 

a  2  i 

u  =  \  f  El  l^4)2  dx  +  i  /  mgU-x)  (|£)2  dx  (8) 

^  o  3x^  ^  o  9X 


t--: 


p  f 


E 

R 

i'.'- 
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A  structural  damping  force  and  a  shear  force  are  assumed  to  be  the 
two  major  nonconservative  forces.  Furthermore  it  is  assumed  that  this 
structural  damping  is  a  viscous  damping  force  of  the  form 


damp 


D 

u  at 


2?1  “1  m  It 


(9) 


From  momentum  considerations,  the  shear  force  at  the  end  of  the  line  may 
be  given  as  * 

R  *  mf  V  [ft  U)  +  V  fx  (10) 

It  then  follows  that 

l 

«W  •  -  /  2?^  m  |£fiy  dx  -  mf  [|£  {l)  +  V  f£  (£)]  6y(JL)  (11) 


*  See  Figure  3. 


To  derive  the  general  equations  of  motion  from  Hamilton's  Principle, 
equations  (7),  (8)  and  (1  )  are  substituted  into  equation  (1)  first.  The 
procedures  of  the  calculus  of  variation  are  then  performed.  Appendix  I 
details  the  mathematics  involved.  After  considerable  algebra,  the  partial 
differential  equation  of  motion  is  found  to  be 

EI  a/"  +  (m  +  mf }  3t^  +  2mf  V  ^  +  mf  V2 

+  2^  m  -  mg  ^  [U-x)  =  0  (12) 

and  the  boundary  conditions  are* 


El  (Z)  =  0 

3x 

(13) 

£  <0)  *  ° 

(14) 

El  ^  (l)  -  0 

3x 

(15) 

o 

n 

o 

(16) 

1 . 2  The  Nonlinear  Tension  Effect** 

A  very  important  physical  characteristic  of  the  tube  undergoing  flutter 
is  its  large  amplitude  periodic  motion.  This  can  only  be  explained  by 
taking  nonlinear  effects  into  account  because  no  linearized  model  can  predict 
such  a  periodic  solution. 

*  In  Part  III  these  are  changed  to  those  appropriate  for  simply  supported  edges 

**  Up  until  now,  it  has  been  believed  that  this  is  the  only  study  on  the  non¬ 
linear  tension  effect  on  the  dynamics  of  a  fluttering  propellant  line  (tube). 


Imagine  that  the  tube  is  deflected  as  shown  in  the  free  body  diagram 
for  a  differential  length,  ds  (Fig.  4).  For  small  displacements,  the  approxi 
nations 

0.|  ^  sin  0.]  'v  9y/3x  (17) 

S2  ^  sin  02  •v*  (§£)  dx  (18) 

may  be  used.  Summing  forces  in  the  y-di recti  on  due  to  tension  yields  a  net 
force  of 

To  I  -  To  <!£  +  J  ■  -To  J  dx  <’9' 

for  a  tube  element  of  differential  .length,  dx.  Assuming  that  the  propellant 
line  is  Hookean  ideal  elastic,  then 

T0  -  axa  (20) 

0Y 

E  =  -r-  (21) 

ex 

where  ex  =  T  (22) 

The  elongated  length  can  be  written  as 


Neglecting  higher  order  terms,  equation  (22)  may  be  written  as 

/  n  (f)2i  ix  - 1  . 

Ex  =  i  21  /  dx 


Combining  equations  (25),  (21)  and  (20)  give 


The  coordinate  system  (C,  <}>)  and  time  t  chosen  to  nondimensionalize  the 


equations  of  motion  are  defined  as 


K  =  x/a 

*  =  y/r,  (r2  »  R2/2)* 


i 


(29) 


t  (30) 

With  the  above  dimensionless  quantities  substituted  into  equation  (28),  the 
equation  of  motion  is  found  to  be 


11 


_  2 

2?, to-,  mi, 

v  s  -  *  (dampinq  parameter) 

vtl(m+mf) 

2 

and  the  fact  that  for  a  thin  cylindrical  tube,  I  =  ar  ,  equation  (30)  can 
be  expressed  as 

^  +  ^l!r  ♦  !M  /  (t)2d£-  'J=0  (3,) 

Equation  (31)  has  associated  boundary  conditions 


at  5  =  0 

* 

o 

II 

(32) 

at  %  a  1 

o 

P  U ) 

3? 

2^-o 

3CJ 

(33) 

2.2  Galerkin's  Method 

The  partial  differential  equation  discussed  in  the  previous  section  may 
be  reduced  to  a  set  of  ordinary  differential  equations  by  using  Galerkin's 
method.  It  is  assumed  that  the  solution  is  of  the  form 

m 

<j>  (?,  t)  *  z  Yi,  (?)  qk  (T)  (34) 

k=l  k 

where  qk  (t)'s  are  unknown  functions  of  time,  and  yk  (?)'s  are  assumed 
functions  satisfying  all  boundary  conditions  of  the  problem.  For  the  uniform 
cantilever  beam  considered  in  this  study,  yk  U)'s  are  taken  to  be 

Yk(0  ■  cosh  XkC  -  cos  Xk5  -  ak[sinh  X^  -  sin  X^]  (35) 

Values  of  Ak  and  ak  are  given  in  Table  1.  Substituting  equation  (34)  into 
equation  (31),  and  using  the  notations 

*  =  3.52  yil/mT 


12 


3qk(T) 
3x  “ 


qk(T) 


The  following  equation  is  obtained 


3YkU> 

3C 


= 


(36) 


Z  qk(x)  Yk1V(5)  +  Z  qk(x)  Yk(S)  +  2 6^1  E  Yk(C)  qkU)  +  6  Z  qk(x)  Yk(0 
=1  K  K  k=l  k=l  K  k=l  K  K 

(37) 

m  i  1  in  i  p  m  h 

+  v  I  Yk(U  qk(T)  -  j  [f  (  Z  qk  YkU))  <£]  *  qk(x)  Yk(5)3  =  0 


The  integral  in  the  last  term  of  equation  (37)  may  be  re-written  as 


1  m 


m  m 


1  . 


/  (  2  qk(T)  YkU))  d?  =  Z  I  [qh(T]  q,(x)  •  f  Yh(C)  yAZ)  d?]  (38) 

o  k=l  K  K  h=l  i=l  n  1  o  "  1 


The  mode  shape  integral  is  evaluated  in  Appendix  II  with  its  values  tabulated 
as  in  Table  2.  Equation  (38)  then  becomes 


1  m  i  n  mm 

f  (  I  qk(T)  Yk(Or  dC  =  Z  2  chi  qh(x)  q.(x) 
o  k=l  K  K  h=l  i=l  ni  n  1 


(39) 


Applying  Galerkin's  method  to  equation  (37)  yields 


m 


Z  ([/  YjYk  1V  dC]  qk  +  [f  YjYk  d£]  qk 


+  28  4T  [f  YjYk  d£]  qk  +  62  [  /  YjYk  d£]  9k 


■  -I  m  m  i 

+  v[/  YjYk  «C]  qk  -  7^  ^  qh(t)  q,(T)  ch1{/  YjYk  df]  qfc} 

(j  —  1 »2, . . . ,  m) 


1 


=  0 


(40) 


To  simplify  equation  (40),  define 


ajk  *  /  d5  (41 1 

bjk  ■  /  t}\  «  <«> 

a^  and  b^  have  been  evaluated  in  the  same  manner  as  c^.  Their  values  are 
also  tabulated  in  Table  2.  Furthermore, 


/  ^k « ■  {jk 


1 


{  '  sjk  \* 


(43) 

(44) 


where  6^  is  the  delta  function.  Equation  (40)  may  be  further  reduced  to 


m 


(2,  {sjk  \  *  I26  *  ajk  +  v6jkJ  qk 


4  p  -.mm 

+  l«4k  xk  +  8  bjk  -  7  ,2,qh  q1  chl  bjk)  qk> 


■  0 


(j  —  1 ,2,...,m) 


(45) 


2.3  Numerical  Methods 

The  time  history  of  the  tube  motion  may  be  determined  by  solving  equation 
(45)  numerically.  For  simplicity,  first  consider  a  two  mode  analysis  for 
equation  (45)  without  the  nonlinear  term,  that  is 


i,  {<5jk  qk  +  t28  ft  ajk  +  6jk  ^  qk  +  ^jkxk  +  8  bjk^  qk}  =  0 


(j  =  1,2) 


In  another  form. 


q-|  +  [20  y/u  a-|i  +  v]  qj  +  [Xj4  +  02  b^J  q1  +  20  47  a12  q2  +  0^  b12  q2  *  0 


20  47  a2-|  q i  +  3  bgi  Q]  *  q2  [20  47  a22  +  v]  q2  +  [X2  +  0  b 22 ]  q2  =  0 


If  one  defines 


9 3  =  q4  =>  ql  “  q3 


q4  =  q2  =>  q2  =  q4 

equations  (47)  and  (48)  may  be  written  in  matrix  form  as 


qg  ~ ( X -j  +0  b-j-j)  -0  b-|2  _(2 047  a-j-j+v)  —  20»^T  a^  qg 


.q4J  L  "^21  -(^2+82b22)  -2041  a21  -(2047  a22+v) 


In  short. 


q  =  A  q 


S(x)  *  aU+Ax)  -  q(T, 
^  —  At 


q(x+Ax)  =  Aq(x)  Ax  +  q(x) 


for  Ax  sufficiently  small. 


The  fact  that  equation  (54)  is  easily  expandable  for  an  m-mode  analysis 


makes  it  the  preferred  form  for  writing  the  equations  of  motion. 

In  the  case  of  four  mode  nonlinear  analysis,  the  nonlinear  term  of 
equation  (45)  alone  generates  256  terms.  However,  only  20  of  them  are  unique. 
The  size  of  the  matrix  can  thus  be  reduced  to  (8,28).* 

The  error  of  the  approximated  solution  can  be  reduced  by  using  the  more 
accurate  finite  difference  relation  ** 


q(x+Ax)  ^  q(x)  +  At(1  +  ^  V  +  5/12  +  _ )  q(x) 

(55) 

where 

Vq(x)  =  q(x)  -  q(x-Ax) 

(56) 

and 

V^1  q(x)  =  Vr  q(x)  -  Vr  q(x-Ax) 

(57) 

A  listing  of  WATFIV  computer  programs  for  four  mode  linear  and  nonlinear 
analysis  is  given  in  Appendix  IV.  They  are  based  on  the  algorithms  discussed 
in  this  sections  and  improved  by  equation  (55).  Initial  conditions  and  data 
used  in  the  program  are  given  in  Table  3. 

The  above  approach  to  determination  of  nonlinear  limit  cycle  oscillations 
is  similar  to  that  employed  by  Dowell  for  plates  and  shells  (7). 


*  Note  that  the  vector  and  the  matrix  have  to  be  conformal. 

**  The  formulas  given  can  be  found  in  Reference  6. 


3.1  Comparison  of  Numerical  and  Previous  Solutions 


The  three  most  interesting  questions  posed  are,  for  a  given  mass  ratio 
and  damping  ratio, 

i)  What  is  the  fluid  velocity  at  which  the  tube  begins  to  flutter? 

ii)  What  is  the  flutter  frequency?  and, 

iii)  What  is  the  limit  cycle  amplitude  of  vibration? 

Answers  to  all  the  questions  above  may  be  obtained  from  the  raw  computer 
data.  The  most  primitive  method  is  the  "binary  search"  method.*  The  computer 
program  is  run  at  several  values  of  the  fluid  velocities  until  the  flutter 
behavior  is  observed.  Binary  search  is  then  applied  to  determine  the  fluid 
velocity  at  which  the  system  is  neutrally  stable. 

Another  method  is  to  process  the  raw  data.  Frequencies  and  damping 
ratios  of  the  different  vibration  modes  are  plotted  against  the  nondimensional 
fluid  velocity  on  the  same  graph.  Flutter  velocity  and  frequency  may  be 
readily  detected  from  the  graph.  Examples  of  two  and  four-mode  linear 
analysis  are  shown  in  Figures  5  and  6. 

To  make  sure  that  the  numerical  results  obtained  reflect  the  actual 
dynamic  instabilities  of  the  tube,  rather  than  numerical  instabilities,  the 
results  from  the  four  mode  linear  analysis  are  double  checked  against  the 
previous  solutions  derived  by  Alan  S.  Greenwald  and  John  Dugundji  (1). 

Although  the  approaches  used  are  significantly  different,  the  two  results  are 
found  to  be  In  very  close  agreement.  They  are  shown  in  Figures  7a  and  7b. 

*  Concept  of  Binary  Search  is  given  in  detail  in  Reference  8. 


17 


3.2  Theoretical  Limit  Cycle  Predictions 

From  the  nonlinear  model,  one  can  predict  the  behavior  of  a  cantilevered 
tube  conveying  fluid  at  low  fluid  velocity.  For  a  given  mass  ratio,  the 
amplitude  of  the  vibrations  is  expected  to  decay  to  zero  with  time.  However, 
above  a  certain  critical  velocity,  the  amplitude  would  increase  first  as  if 
there  is  no  upper  bound.  After  a  few  cycles,  the  amplitude  of  the  vibrations 
drops  back  by  a  fraction  and  stays  at  a  constant  magnitude  indefinitely.  This 
is  called  a  limit  cycle.  The  amplitude  of  the  limit  cycle  increases  with  the 
fluid  velocity.  Graphical  presentations  of  the  relationship  between  the 
amplitude  and  the  fluid  velocity  are  given  in  Figure  8  and  9  for  mass  ratios 
of  0.1  and  0.5  respectively. 

First  and  second  mode  frequencies  are  identical  at  the  threshold  of 
instability.  Beyond  that,  both  frequencies  increase  with  the  fluid  velocity. 
This  is  shown  in  figures  10  and  11  for  mass  ratios  of  0.1  and  0.5  respectively 

3.3  Mode  Shapes  and  Stress  Distributions 

Mode  shapes  of  the  fluttering  tube  can  be  obtained  by  holding  t  in 
equation  (34)  constant.  Altogether  five  values  of  t  are  judiciously  chosen  to 
represent  equal  time  increments  over  one-half  of  the  period  of  the  limit  cycle 
oscillations.  The  mode  shapes  of  the  tube  at  these  instances  are  plotted  on 
the  same  scale  for  ease  of  comparison.  Figures  12,  13,  14  show  mode  shapes 
of  tube  of  mass  ratio  0.1  for  three  different  fluid  velocity,  B.  Figures  15, 
16,  17  show  the  cases  where  mass  ratio  is  0.5,  and  the  mode  shapes  are  at 
1/2,  5/8,  3/4,  7/8  and  the  end  of  the  period. 

Stress  distribution  may  be  determined  by  using  the  following  equation 


or  ax  =  ~1 - 1 

x  r  k 


Er 


m 


=  7T  kf1  qk(T)  Yk(c) 

-2  m 


(59) 


Er  111  7 

=  ~y  £  qk(x)  •  [cosh  Ak?  +  cos  AkS  -  ak  (sinh  X^E,  +  sin  \kE)] 

&  k-1 


Figure  18  shows  the  stress  distribution  in  the  tube  for  y  -  0.5  and 
B  =  10.4  at  five  different  time  steps.  Unlike  what  was  expected,  the  maximum 
stress  does  not  occur  at  the  root  of  the  tube  (the  clamped  end).  It  appears 
to  be  shifting  along  the  length  of  the  tube  during  the  entire  cycle  of 
oscillation.  This  is  because  the  maximum  curvature  of  the  tube  does  not 
necessarily  occur  at  the  root  of  the  tube. 


*  o  is  in  psi  in  Fig.  18. 


Section  1 


Equations  of  Motion 


The  partial  differential  equation  of  motion  of  a  simply  supported  tube 
is  the  same  as  the  one  for  a  cantilever  tube,  i.e.,  equation  37.  The  only 
difference  is  the  set  of  boundary  conditions 


at  £  =  0  >  $  =  0 

4=o 


(60) 


at  £  =  1,  <p  =  0  'j 

4=o 

It  is  again  assumed  that  the  solution  is  of  the  form 


(61) 


t) 


m 

I 

n=l 


Y„(S)  q„(T) 


162) 


To  satisfy  the  boundary  conditions,  the  yp(£)'s  are  taken  to  be 

YnU)  =  sin  mr£  (63) 


n  -  1,  •••  j  m 

Since  everything  else  is  the  same  as  for  the  cantilevered  tube  except  for 
Yn(S)»  equation  (40)  can  be  reused.  The  only  things  that  have  to  be  rederived 
are  the  mode  shape  integrals.  Mode  shape  integrals  for  a  tube  simply- 
supported  at  both  ends  are  evaluated  in  Appendix  V. 

As  before,  define 


\jk 


1 

=  / 
o 


(64) 
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bjk  ■  /  « 


i  .  . 
Ch1  =  !  Wi 


and  the  equations  of  motion  become 


m  i  »  > 

{I  6jk  qk  +  [2B  ^  ajk  +  2  v  6jk]  qk 

k4  4  2  i  m  m 

+  [nr  5 ,k  +  b.k  .  J.  ^  ^  %  q.  Ch.  bjfc]  qk>  -  0 


(j  ~  1 »2,  •••  j  m) 

The  computer  program  written  for  evaluating  the  time  histories  of  a  canti 
levered  tube  may  be  readily  modified  to  solve  equation  (67)  by  the  following 


changes 


j  j=k 

0  j^k 


Xk  =  kir  (69) 

A  listing  of  the  modified  program  in  FORTRAN  is  given  in  Appendix  VI.  Other 
alternations  were  made  to  conserve  core  storage  space  and  computer  time. 


Secti 


ion  2 

Predicted  Characteristics 


For  0  £  0  <  it,  the  tube  undergoes  damped  stable  vibrations.  Examples 
of  these  damped  vibrations  are  shown  in  Figures  19  and  20  for  the  two  limiting 
cases.  The  magnitude  of  the  true  damping  ratio,  ej,  increases  with  flow 
velocity.  This  relation  is  shown  in  Figure  21a.  First  mode  frequency,  u-j , 
as  a  function  of  0  is  shown  in  Figure  21b.  The  frequency  of  the  first  mode 
vibration  decreases  as  0  is  increased  n.  In  the  neighborhood  of  0  =  it,  the 
first  mode  damping  ratio  becomes  very  large  and  the  frequency  approaches  zero. 
In  Figure  22,  the  product  is  shown.  Note  that  it  varies  smoothly  with 
0  as  0  passes  through  tt. 

Beyond  0  =  tt,  the  tube  was  found  to  undergo  divergence  by  previous 
investigators.  Ref.  5.  This  is  only  true  for  0  less  than  2ir.  For  tt  <  0  <  2tt, 
the  magnitude  of  q-|(x)  grows  exponentially  with  time  initially,  q^ (t )  then 
oscillates  near  some  constant  value,  say  q*.  These  transient  oscillations 
decay  with  time  and  eventually  die  out,  and  (x)  converges  to  the  steady 
state  amplitude,  i.e.,  q*.  Higher  mode  oscillations  are  not  very  significant 
and  the  amplitudes  decay  to  zero.  Examples  of  such  behavior  for  four  dif¬ 
ferent  values  of  0  are  given  in  Figures  23,  24,  25  and  26.  To  show  the  re¬ 
lationship  between  0  and  the  amplitude  of  q.|(x)  more  clearly,  the  non- 
oscillatory  average  of  q-|(x)  is  plotted  against  time  for  the  four  0's  in 
Figure  27. 

Near  0=  2tt,  all  four  modes  become  dynamically  unstable.  However,  the 
first  and  second  modes  are  dominant.  Four  of  the  relevant  time  histories  are 
illustrated  in  Figures  28  to  31. 
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It  was  speculated  that  if  at  6  =  tt,  the  first  mode  becomes  statically 
unstable,  and  at  8  =  2ir,  the  second  mode  becomes  important  and  changes  the 
static  instability  to  a  dynamic  Instability,  then  at  8  =  3ir,  the  third  mode 
may  become  very  significant  and  drastically  alter  the  time  history  of 
motion  again.  However,  this  is  not  true.  The  third  and  fourth  modes  are 
dynamically  unstable,  but  they  do  not  dominate  the  total  motion.  For  B  >  3tt, 
no  drastic  change  in  time  history  is  found.  It  appears  to  have  the  same 
"pulse"  or  "beat"  characteristics  as  for  B  <  3ir.  These  are  shown  in  Figures 
32  and  33. 

A  stress  analysis  was  done  for  B  =  5.92.  As  was  expected,  the  maximum 
stress  is  in  the  middle  between  the  two  ends.  Figure  34  shows  the  stress 
distribution  in  the  tube. 

Although  there  is  no  drastic  change  in  the  time  history  for  B  >  3tt,  the 
third  mode  does  appear  in  the  deformation  shape  during  the  cycle.  At  8  = 
10.0,  the  deformation  shape  of  the  tube  changes  from  a  two  mode  sine  wave  to 
a  three  mode  sine  wave  and  back  to  a  two  mode  sine  wave  over  a  complete  cycle 
of  motion.  Deformation  shapes  at  20  equal  time  increments  over  the  period  of 
the  modulated  sine  wave  are  shown  in  Figures  35a-35e.  Graphs  showing  stress 
distribution  in  the  tube  for  the  20  corresponding  deformation  shapes  are 
given  in  Figures  36a-36e. 

The  effects  of  fluid  velocity  B  on  the  first  and  second  mode  peak  limit 
cycle  amplitudes  are  shown  in  Figure  37a  and  37b  respectively.  These  further 
illustrate  the  dramatic  changes  in  modal  content  of  the  motion  which  occur 
near  b  =  2tt. 

The  first  correct  linear  stability  analysis  of  a  simply  supported-simply 
supported  pipe  was  given  by  Housner  (9).  The  present  discussion  extends  his 
results  to  the  nonlinear  regime. 


This  study  has  developed  a  mathematical  model  capable  of  describing  the 
large  amplitude  periodic  behavior  of  a  cantilever  or  simply-supported  pro¬ 
pellant  line  (tube)  conveying  high  velocity  liquid.  The  study  has  also 
illustrated  the  importance  of  the  nonlinear  tension  effect  induced  by  bending 
of  the  tube.  Thresholds  of  instability  predicted  by  the  theory  are  in  good 
agreement  with  those  previously  obtained  by  other  investigators.  Plans  have 
been  made  for  setting  up  experiments  to  confirm  the  predicted  magnitudes  of 
the  limit  cycles. 

The  theoretical  model  might  be  improved  by  considering  other  non- 
linearities,  e.g.,  nonlinear  modulus  of  elasticity  or  nonlinear  curvature. 
Instead  of  viscous  damping,  one  might  consider  viscoelastic  damping,  or 
better  still,  a  combination  of  viscous  and  viscous-elastic  damping.  The 
uniform  beam  assumption  could  be  replaced  by  the  thin  cylindrical  shell 
approximation.  For  cases  in  which  the  hanging  tubes  are  very  flexible,  the 
gravity  term  should  not  be  neglected. 
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Table  3 

INITIAL  CONDITIONS  AND  DATA  USED 
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Figure  3.  Tube  with  Flowing  Fluid 
(From  Reference  1) 
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Figure  6a.  The  Effects  of  Velocity  on  Frequency  (u  =  0.1) 

(Two  Modes) 
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Figure  11.  Variation  of  Frequency  Due  to  Velocity 
(y  *  0.5,  Four  Modes) 
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Figure  14.  Tube  Deformation  Shapes  at  Various 
Times  During  Limit  Cycle 
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Figure  17.  Tube  Deformation  Shapes  at  Various 
Times  During  Limit  Cycle 
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Figure  19b.  Second  Mode  Time  History  of  Tube  Motion 

(Damped  Vibrations) 
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Figure  35e.  Tube  Deforamtion  Shapes  at  Various 


Time  Intervals  (8  =  10.0) 
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APPENDIX  I 


Solving  the  Energy  Equation 

Substituting  equations  (7),  (8),  and  (11)  into  equation  (1)  gives 


1 « /  !  ‘"’(f)2  +  +  2V  If  v2  (f>2  +  v2i 

t,  o 


El  +  mg(i-x)  (|^7)2)  dxdt  -  f  f  2?i^m  |^-  3y 

3x  t-|  o  ' 


+  /  mf  [|£  (£)  +  V  (£)]  •  dt  6y  (£)  =  0  (A-l) 


In  the  calculus  of  variations,  <$  may  be  considered  as  a  linear  differen 
tial  operator.  Examples: 
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Similarly,  performing  variation  on  equation  (A-l)  and  then  rearranging  the 
terms  yields 
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t 

1  (A-5 

Due  to  the  tedious  mathematics  involved,  only  representative  terms  of 
equation  (A-5)  will  be  worked  out  completely  and  shown  below. 

Integration  by  parts: 
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(See  Ref.  6} 


<5y  dx 
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To  satisfy  Eq.  (A-5),  the  integrands  have  to  be  zero.  The  sum  of  all  the 
integrands  of  the  double  integrals  give  the  partial  differential  equation  for 
the  problem. 


El  H  +  lm+mf)  Z$  +  2mfV  Mx  +  mfv2  H  +  2nVl  It  "  mg  It  [(*"X'  lf]  =  0 
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Also  the  several  end  point  terms  must  be  zero.  Thus 


(A-13) 


For  a  cantilevered  tube  (A-13)  is  satisfied  by 


dX  >X=2. 


Fcr  d  simply-supported  tube. 
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The  other  end  point  conditions  are 


[EI  ^  -  mg(Jl-x)  |£]  6y 
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For  a  cantilevered  tube,  (A-14)  is  satisfied  by 

Wv.n  *  0 


[EI  ^4  -  mg(i-x)  |£|  -  0 
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For  a  simply-supported  tube. 


1  x=0 ,  l 


(A-14 


It  should  be  noted  that  in  obtaining  (A-14)  from  (A-5)  to  (A-ll),  there  has 
been  a  cancellation  of  terms  involving  the  last  term  in  (A-5)  and  the  first 
terms  on  the  right  hand  side  of  (A-7)  and  (A-10)  for  x-l  to  give  the  same 
form  in  (A-14)  for  the  boundary  conditions  at  both  x-0  and  x=JL 
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Integrals  containing  the  characteristic  functions  and  their  derivatives 
are  referred  to  as  mode  shape  integrals  in  this  paper.  They  can  be  evaluated 
by  the  method  of  partial  integration  along  with  their  orthogonality  proper¬ 
ties.  Special  formulas  for  those  mode  shape  integrals  which  occur  frequently 
in  engineering  application  may  be  found  in  references  10  and  11.  They  are 
given  in  terms  of  an  and  Bn  for  five  common  boundary  conditions:  clamped-clamped, 
clamped-free,  clamped-sup:  jrted,  free-free,  and  free-supported. 

As  an  illustration,  consider 


1  , 

chi  *  f  Yh  <0  Y1  (C)  d£ 


(A— 1 5) 


for  a  cantilevered  uniform  beam.  Combining  formulas  (6)  and  (16)  of  reference 
(7)  yields. 
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8h  -  «i 


The  corresponding  a's  and  B's  can  be  found  in  reference  11.  They  are 


'•TV 


7.8546 

10.9955 
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G7  61536.SILFl!943I'.  =  .r.  82 


♦  + 

♦  NUMERICAL  SOLUTION  FCr!  STATIC  AND  + 

♦  DYNAMIC  INSTABILITIES  OF  A  PROPELLANT  ♦ 

+  LINE.  (FOUR  .NODES  ANALYSIS)  ♦ 

+  NP:  WITH  FORMULA  CF  OPEN"  TYPES  + 

♦  + 


INTEGER 

II 

,JU,:n, 

iJ2, 11 

,12 

FEAI 

CC  ( 

Q 

°  # 

1500) 

F.EA1 

ELMFY 

,hat:i 

S  £  A I 

3  KU 

REAL 

E£T  A 

,ao,ao. 

x  i  E 

P  E  A I 

ci. 

02 

,Q3,Q4, 

11  ,L2 

,vr 

REAL 

05, 

06 

,07, 08 

E  E  A  I 

L3, 

1.4 

E  E  A I 

L  ,Z 

1, 

I ,  TMAX  , 

DFLT , 

All 

FEAL 

All 

* 

i  ^ 

22 , B 1 2 , 

522,221, 

?  E  A 1 

A  13 

* 

f  *• 

1 4 ,  .■»  2  - , 

A24 

REA  I 

A  31 

t  * 

?  2  ,  A  3  2  , 

A  34 

FEAL 

A  41 

$  *- 

4  2  ,  A  4  3 , 

A44 

FEAL 

E  13 

9  - 

1 4 , 3  <4  3  , 

a  24 

REAL  331,522,933,534 

SEAL  E  4 1 ,542,943,944 

F.EAL  K  (8,8)  ,  A  2  (3,1500)  ,AI  (1500)  ,V  (8) 

HE  A I  TQ 1  (5C0)  ,ri2  (5GC)  ,11’  ( 5 C 0 ) 

SEAL  1C3  (500)  ,  IC4  (5C3) 


Q1  ,  02,  Q5,  ...  ARE  INITIAL  CONDITIONS 

CO  (\  )  •  S  Ar£  PSEUDO  INITIAL  CONDITIONS 

11,  L2  k 3. £  ”  DDF  SHAFT  PATAMEIEnS 

L  *  LENGTH  OF  PIPE 

Z1  =  CRITICAL  DAMPING  RATIO 

f  *  HASS  CF  PIPE 

HP  =  HASS  OF  FLUID 

VF  =  FLUTE  VELOCITY 

I  =  I N  F  R  T I A 

E  =  MODULUS  OF  FLASTICI",Y 

All,  A  12,  . ..,  B11,  512,  ...ARE  CODE  SHAPE  INTFGP.ALS 


REAL  (5, 10C0)  01,02,03,04,05,06 
REAL  (5,1000)  07, 03, LI ,L2 ,L3,L« 

FC2KA1  (6F1T.4) 

F.EAC  (5,1CCi)  A11.A12,  A?  1  ,  A  22 ,  BIT  ,  C 1 2 , 32 1 , 32"> 
FEU  (5  ,  1  0  C  1 )  A  13,. A  14,  A  23,  A  24, 513, £14,323,324 
REAL  ( 5  , 1 0  C  1 )  »21,A3?,A33,A34,141  ,i42,A43,.A44 
PEAL  (5 , 1 r  Cl }  531,332,233,534,541  ,EU2, 343, 344 

FC.SRU  (H  FI  0.4) 

FEU  ( 5 , 1  1  11,  I  "I  C  =  3 f|  q )  VF,  1,71,1,  CELT, THAI 
I  FOM"”  (F12.4,5E10.U^ 
hE.AE  (5,  KC-)  f , I , I 

FCPMAT  (4 [ 10.3) 


BET  1  =  Dlf.TNSISNLIFS  FLUID  1 
HU  =  .HASS  PATIO 


c 

c 


MJ 


105 


=  VISCCUS  DAF.PI3G  COEFFICIENT 


32  BITS  =  V?*JQ.tT  (:i?*L**2  ,/(E*I)) 

33  BC  =  SF/dr  +  Y) 

34  NU  =  7.0u*Zl*SQrT(1.-^U) 

C 

35  WHI1F  (6,1012) 

3*  WHITE  (6 ,  1  CO  3)  ?.  1 1 ,  A 1 2 ,  .A  1  3 ,  A  14  ,  A2  1 ,  A22 ,  A23 ,  A24 

37  1003  FCSKAT  (•  • , • ?  (1 , J)  =  ‘ , 4F2C. 4,/, •  •  ,  •  A  (2,  J)  =  •  ,  4F20 .  U) 

33  WFITE  (6 , 3  C  1 1 )  f  3 1 ,  A32 ,  AJ3,h3U  ,? 4  1 ,  A42 ,  A43 ,  A44 

39  3011  FCSEAT  (•  • , • A  ( 3, J) = • , 4F2 0 . 4 ,/ , •  •  ,  •  A  (4 , J)  =  • , 4320 . 4) 

40  WHITE  <6,3C01)  E11f312,E13,B14,B21  ,322, 323, P 24 

41  3C01  FCEKAT  (•  •  ,  •  B  (1 ,  J)  =  •  ,4F2C.  4,/,  •  •  ,  •  3  (2,  J)  =  '  ,  4?20. 4) 

42  WHITE  (6,3111)  E 3 1 , B 3 2 , E3 2 , P 34 , 5 4 1  ,  P42 , 34 3 , Bu4 

43  311  1  FCFKJ.7  ('  '  ,  •  3  (3  ,  J)  =  •  ,  4F2C.  4  ,/,  •  »  ,  *  3  (4  ,  J)  =  • ,  4F20.  4) 

WHITE  (6,1004)  Cl , 02 , Q3 , 0 4 , Q5 , Q6 , 07 , 08 

45  1  CO  4  JOSKJT  (.  •,*0(K)=',4F20.4,/,f,X,4F2C.4) 

46  WHITE  (6 , 1  C05)  BETS, ‘'0, NO 

47  10C5  FCSM?"  (•  •,»3ETA*l,1?15.4,5T,*tCU»f,1F15.4,5X,»S'J»'f1Fl5.4) 

48  WHITE  (6,1006)  VF,L1,L?,L2,L4,I,T1 

49  1006  FCSEM  ('  ' ,  •  VF=*  ,  1F1  5.4,  £X,  •  LI  ...  ASS' ,  4  f  10 . 5  , /,  •  ','L=', 

15X,  '  21  =  ’  ,  1 F 1 0 . 4) 

5C  W-I1E  (6,iC07)  M,1F,I,2 

51  1007  FCSKAr  (•  ','E,LF,r,ii  SEE 4620. 3) 

52  WSITE  (6 , 1  C06) 

53  1008  fC3K.iT  (//,30X,  'MAT BIX  IS',/) 

C 

54  EMU  =  £QFI  (MU) 

C 


c 


55 

DC  (1,1) 

=  C.O 

56 

DC  (1,2) 

=  0.0 

57 

DC  (2,1) 

=  0.0 

58 

DC  (2,2) 

=  0.0 

59 

DC  (3,1) 

=  C.O 

60 

DC  12,2) 

=  C.O 

61 

DC  (4,1) 

=  C.O 

62 

DC  (4,2) 

=  C.O 

63 

DC  (5  ,1) 

=  C.O 

64 

DC  (5,2) 

=  C.O 

65 

DC  (6,1) 

=  0.0 

66 

DC  (6,2) 

=  V  .  0 

67 

DC  (7  ,1) 

=  C.O 

68 

DC  (7,2) 

=  C.O 

65 

DC  (6,1) 

=  C.O 

7C 

DC  (8,2) 

=  C.O 

71 

K  (1,1)  = 

C.O 

72 

*(1,2)  = 

o.c 

73 

Ml, 2)  = 

C.O 

74 

*(1,4)  = 

0.0 

75 

*(1,6)  = 

1.0 

76 

*(1,6)  = 

0.0 

77 

*(1,7)  = 

C.O 

78 

*  (1,6)  = 

0.0 

79 

*(2,1)  = 

C.O 

in  in  n  (n  c'i  r»>  «'i  to  m  i»)  m  m  rr  rr  cr  -C7  a-  a-  ^  a  1/1  in  uun  ut  mm  u)  vcvo  'o  i>r>  r*-  r~  r~  r*-  p*  a.>  a>  a>  a’ 


z«) 

.*) 

,D 

.*) 

z1) 

,2) 

#2) 

z^) 

z5) 

zE) 

z'/) 

zf) 

zl) 
z2) 
z  2) 
z<0 

z  2) 

,t) 

z'/) 
zE) 
zl) 
z  2) 
z  2) 
,“) 
z  2) 

z  E) 

z7) 

zf) 

zD 

z2) 
z  2) 
,L‘) 
z  2) 
z  E) 

z'O 

zE) 
zl) 
z2) 
z  2) 
z«> 
z  2) 

z  E) 
Z  V 

zf) 

zD 

z2) 
z  2) 
z«) 
?,£) 
9  z  E) 
5,7) 
9,E) 


(L1**4.*-33TA**2.*3l1) 

E£?A**2 . *31 2 

BETA  *  *  2 . *313 

B37A**2.*S14 

(2.*JU?A*H:iH*A1  1  +  NO) 

2.  *3ETA*nnt’*A12 

2. *3ETA*3«0*A13 

2.*3EIA*3f!U*Al4 

EETA**2.*B21 

(L?**4.*SE?A**2.*322) 

££TA**2.*B23 

E2IA**2.*B24 

2 .  *5  21?.* -.*’’*  A21 

(2.  *BETA *?VJ*A22*NU) 

2.*3ETA*?.?TJ*A23 

2.*BEIA*3HU*A24 

EETA**2 .*331 

EEIA**2.  *B.?2 

(L3**4.*3E!A**2 .*333) 

EE T?  **2 . *3  34 

2.  *3E1 A*?.  MU*  A3 1 

2. *  SET  **£*3* A3 2 

(2.*BET^*EM.  U*A33  +  NU) 

2.  *  3  £  I A  *  1 !’  J  *  A  3  4 

EErA*»2.*B4 1 

B2TA**2. *342 

ESTi**2.*B43 

(L4*  *4.  +  BETA**  2.  *344) 

2.  *c£T»*2MrJ*A41 

/  .  *  3  ET  A  *  EMU*  A4  2 

2  .  *  3  £T  A  *  2  M  rJ*  A  4  3 

(2.  *3ETA*E  iU*A4  4+UU) 


1ZI  (5 , 1  C09)  ((K(II,JJ),  J  J=  1 , 6 ) 

rtf' AT  C  ',3? 12.4) 


133 

AC(1,M)  =  01 

139 

AC(2,>M)  =  Q2 

I4C 

AC(3,N1)  =  03 

141 

AC  I  4  ,  N  1 )  =  G 4 

142 

AC  ( 5  / Ti  1 )  =  Q5 

143 

AC  » €  #  12 1>  =  06 

144 

ACP,M)  =  Q7 

145 

AC  (6#N1)  =  C3 

146 

ICC 

C 

CONTINUE 

147 

DC  1C1  11=1, 

146 

V(I1)  =  AC  (1 1 

149 

101 

c 

CONTINUE 

I5C 

AT  ( K  1)  =  T 

151 

NNN  =  N 1  +  2 

152 

HI  =  N 1*1 

MATXIX  MULTI PL  ICAIICN 


5  3 

DC  1 C  2  13=1,8 

54 

DUMMY  =  0. 

55 

DC  1C?  12=1,9 

56 

MATH  =  X  (I  3,12)  * V  (12) 

57 

DUMMY  =  DUMMY  +  MAT  M 

58 

102 

CONTINUE 

59 

DC  ( I 2 , S N N )  =  DUMMY 

6C 

AC  C 2, HI)  =  CELT*  (1 ,9167*LQ (I3,NNN) 
2  +  0.4167*00  (13,  NSN-2)  )  ♦  V  (13) 

61 

103 

r 

CONTINUE 

62 

L 

I  =  T+DELT 

63 

IF  CMAX.GI.T)  GO  TO  100 

64 

KETIE  (6 , 1  C  10) 

65 

101C 

FCFKAI  (//,5X, »T,01,Q2,Q3,Q4  AEE  AS 

1011 

C 

C 

C 

C 

C 

C 

C 

C 

c 


»CT£:  ONLY  EVSEY  FIFTH  VALUE  15  FEINTED 

NEITP  |6 , 1  Cl  1 )  (AT(J)  ,  (AO  <IZ,J)  ,IZ=1,4)  ,J=1,N1,5) 
FC5KJI  ('  «,5i15.U) 


FLCTTISG  FUNCTION 


NCIE:  ONLY  EVSEY  FIFTH  VALUE  IS  PLOTTED 


68 

t;  3  =  .c 

65 

DC  1  C  4  N  2  =  1 , V 1 

7  C 

S3  =  N3  +  1 

71 

1T(N3)  =  AT  (M2) 

72 

TQ1  (N3)=A0  ( 1 , 13 2 

73 

?C2  <*?)  =  AO  (2, 

74 

TCiC?)  =  AC(3, 

101  3 


1 0 1  a 


1015 


1016 


:c«  <F3)  =  AO  (i»,N2) 
CCJSIVCE 

KHIE  (ft, 1  C12) 

FC’f*;T  C*  1  •) 

CALI  W  FLOI  1  (TT,TQ1,N3 

asm  (6,1013)  T”  ( 1 ) , 

FC5R/I  (40a,  •  DISPLACES 
1  *  T  =  •  ,1F7.4,10X,  'T*AX=' 
W?Ilf  (6,7012) 

C.UI  KFLPil  (TT,Tg?,N3 
WpIII  .  <6,1C14)  Hp), 
FCASAT  (40X,  •  DISPLACE 
1  •  T=  •  ,1f7.4,19X,  '?U.\X=' 

asm  (ft,ici2) 


CALL 

*  FLO 

11 

(i: 

,IC3,N3 

v»:-. 

TIE 

(6, 

1015) 

77(1)  , 

FC 

5KA. 

T  (4 

cx. 

1  r  r 

splace;: 

•I 

-  i 

, 

1F7. 

4 , 1 0  X , 

•ri?.x=' 

W  F 

III 

(5, 

1C1 

2) 

CA 

LI 

m  ELC 

(IT 

,mOU,N3 

Wr 

III 

(6, 

1  C 1 6) 

'"1(1)  , 

FC3KA 

-  (4 

A  «F 

'  DI 

SPLACE;': 

' : 

s  1 

, 

1*7. 

4,1  OX, 

• :nax=' 

he 

(*/ 

1C1 

*) 

GC 

1C 

99ft 

CC 

Nil 

»«  f.  7 

,12, 'DISPLACE WENT') 

T.1A  X,  DELI 

ENI  ( 1 )  AS  A  FUNCTION  OF  ?AU',//,20X 
, 1 ?7.  4 , 1 0X, ' CELT= '  ,  1 ?7. 4) 

,  12, 'displacement') 

IS AX, DELI 

ENT  (2)  AS  A  FUNCTION  OF  TAU',//,20X 
,  1F7.4,1CX,' D£L?= ' , 1F7.4) 

,  1 2 , »  etsflacssent') 

TKA  X, DELT 

ENI  (3)  AS  A  FUNCTION  C?  ?AU',//,20X 
, 1 F7. 4 , 1  OX , ' DDLT=  * , 1F7.4) 

, 12 , »  DISPLACEMENT ' ) 

ISA  X , DELI 

EMI  (4)  AS  A  FUNCTION  OF  T\U',//,20X 
,i?7.4,10X,» DELI-' , 1F7.4) 


SEMFY 


076153b. SZLFFD  ?IX=A32 


♦  NUMERICAL  SOLUTION  FCr  STATIC  AND  + 

♦  EYK&2XC  INSTABILITIES  OF  A  PR Jr ELL A NT  + 

+  LINE.  (FOUR  MODS S  ANALYSIS)  «• 

+  N3:  WITH  FJRMMLa  CF  OPEN  TYPES  + 

♦  NONLINEAR  AND  IMPROVED  + 


INTEGER  IT 


INT5GE 

.  F  I A  , 

Id  , 

ID1, 

T 

X 

D2 

INTEGER  11, 

12, 

13, 

14, 

15 

,  16 

INTEGE 

S  J I 

INIIGEF  J1, 

J2, 

03, 

J4, 

J5 

,  Jfc 

PEAL 

DUMMY, 

M  ATM 

F,  £  A1 

EMU 

REA  I 

EE I A  ,  MU,  HU, 

M  ,  KF 

t 

I,  E 

FI  A I 

Cl,  C2, 

Q3, 

Q4 

•  Q5 , 

06,  Q7,  08 

F.EAI 

LA M  (4)  , 

VF 

.  ^ 

S  E  A I 

It  .21, 

"n  m 

-  F  - 

MAX 

,  DSL 

— 

K  i  .»  1 

DEL  (4,4) 

PEM 

CO (8,3200) 

REM 

A  (^  r  «)  , 

3  (4 

,4) 

,  C(4 

#4 ) 

F.EM 

K (8, 28) 

,  AC  (3, 

3200) 

f 

AT  (3200),  V  (?8) 

REAL 

TCI  (500)  ,  r 

C2  (*C0)  , 

T 

Q3  (50C)  ,  TQ*  (500) 

PEM 

IT  (500) 

Cl,  C2, 

Q  3, 

•  • 

.  A  *4 

il 

INITIAL  CONDI'- 10 US 

EO  CO  'S 

API 

?s 

It'  DC 

IK 

ITIAL  CONDITIONS 

1AM  (K) 

?.  RE 

MOD 

1  SHA 

t  L 

PAPAL  El  EPS 

VF  =  FL 

UIC 

VSL 

OCITY 

L  =  LEN 

GTH 

0? 

PIPE 

Z1  =  CRITIC 

*  L 

DA-MPI 

NG 

RATIO 

K  =  MASS  CF  PIPE 
f?F  =  MASS  OF  FLUID 
I  =  INERTIA 

E  =  MODULUS  ')?  ELAEIICI-’Y 

A  (K)  ,  3  ( N )  ,  C  (U)  ARE  MODE  SHAPE  INTEGRALS 


REAL  (5,1  COO) 

FC3KSI  (8F10.4) 

REAL  (5, 1005) 

REAL  ( 5, 1 C  C5) 

HEAL  ( 5 , 1 C  C5) 

PEAL  (5,1 005) 

FCR  EAT  (4F1Q.U) 

Ff'iL  ( 5,  1C  10) 

EC  PE  SI  («E 10.3) 

REAL  (5, 1 :)  1 1 ,  r'JC=  ■)  R9) 
FORK  AT  (F12.4,EF10.4) 


Cl,Q2,03,Q4,c5,06,Q7,0S 


(LA  I  (II),  I1  =  D4) 

(  (A.  (II  ,J1)  ,  J1=1,L)  ,  11*1,0) 
((3(11, JD  ,  J1=1,4)  ,  11*1,4) 
((C(I1  ,J1)  ,  J1=1,4)  ,  11*1,4) 

f,:if,i,e 

VF,I,Z1,T,CELT,TMAX 


=  0 1 UNLESS  FLUID  VELOCITY 


F£I  n 


29 

30 

31 


L* 


BETA  =  V?*SQHT (M?*L**2./ (E*I) ) 
MU  =  ilt/(ZS  +  ») 

NU  =  7 . 0***Z1+SQF T  ( 1 . -MU) 


32 

C 

W  7  IT  I  (6/1015) 

33 

1015 

FCSM  AT  ( *  1  • ) 

34 

WSITE  (5/TC20)  ((.\( 

35 

ro2C 

F07"A1  (•  '  /  *  A  (I ,  J)  = 

36 

WHITE  (6,1025)  (  (3  ( 

37 

1025 

FCHMA7  (•  ,,*6(I,J)  = 

33 

WHITE  (6,1030)  (  (C  ( 

39 

1 03C 

F0SM*7  (•  •/•C(I/J)  = 

4C 

WHITE  (6,1035)  gl/0 

41 

1035 

FO  5  M  AT  (•  '  ,  '  Q  ( t. )  = '  , 

42 

WFIT2  (6,1040)  6STA 

43 

1 04  C 

FOFMPT  (•  '  , 1 3£TA=' , 

44 

WHITE  (6,1045)  VF,L 

45 

104  5 

FOHMAT  ('  •,'7F=,/1? 

46 

<  *  L  A  fc  EA  (N)  =  '  ,  4?1  C  .  u) 
WHTII  (6,1050)  H ,  K  F 

47 

1 05  C 

FC3KAT  ('  ','K,KF,I, 

43 

W“TIx  (6,1055) 

49 

1055 

FCSM 'I  (// /3CX,  •  MATS: 

5C 

WHITE  (6/1C56) 

51 

1056 

FCSM  AT  ( '  • , '  (DUE  TO 

52 

C 

*  INSTEAD  Cf  (3,28) • , 

EMU  =  SC HI  (MU) 

C 

c 

c 

c 

c 

c 

c 

c 

c 


♦ 

♦ 


3 MTS I S3  OF  MAT SIX 


♦ 

♦ 


NCTl:  IELTA  FUNCTION 


53 

EC  11C  I D 1= 1 , 4 

54 

EC  ICO  ID2=1 , 4 

55 

EEL  (ID1,ID2)  =0.0 

56 

1  co¬ 

CC  ill  INEZ 

57 

DEI  ( ID  1 , ID 1)  =  1.0 

53 

lic 

CONTINUE 

c 

59 

EC  120  14=1,8 

6C 

DC  (14,1)  =0.0 

61 

DC  (14,2)  =0.0 

62 

120 

CONTINUE 

C 

C 

KCIi:  EO-LOCP 

C 

63 

IC  =  4 

64 

DC  1  4  C  I A  =  1 , 4 

65 

EC  120  13=1,23 

66 

K  (15, IE)  =  0.0 

67 

13C 

CONTINUE 

68 

IC  *  IC  ♦  1 

69 

K(IA, TC)  =  1.0 

7C 

14C 

CONTINUE 

C 

■A 


7 


■Psf 

/y 


* 

jf  .#■ 

yj? 

// 

•r 


"5  s>~ 


KC1 1 : 


ENI2IF5  WITH 


COGNIZED  EAITERU 


DC  160  15=5,8 

IT  =  15-4 
DC  150  05=1,4 

K(T£,05)  =  -(DIL(IT,05)*LAZ(IT)**4.043E7A**2.0*B(IT,05)) 

CONTINUE 

DC  155  03=5,8 

07=  J3-4 

K  (15,03)  =  -  (2. 0*32TA*3MU*A  (IT,OT)  ♦DOL  (IT,OT)  *NU) 
CCKIIHUE 


K 

K 

K 

K 

K 

K 

K 

#E 

K 

K 

*E 

K 

HE 

K 

K 

K 

K 

K 

K 

*3 

K 

K 

K 

K 


15,9) 
15, 1C) 
15,11) 
,15,12) 
[75,13) 

15. 14) 

15.15) 

1 1.3) * 

15.16) 

15.17) 

12.4)  * 

75.18) 

17.4) * 

15.19) 
15, 2C) 

15.21) 

15.22) 
-5,2  3) 
-5,24) 

71 .4)  * 
15,25) 
35  ,26) 

15.27) 

15.28) 


=  B  (IT, 
=  E  (IT 
=  E  (IT 
=  E  (IT 
=  BC" 
=  3  (27 
=  5  (IT 
(C  (1,2) 
=  3  (IT 
=  3  (IT 
(C  (1,2) 
=  3  (IT 
(C  (1,3) 
-  D  (IT 
=  3  nr 

=  °  (IT 
=  E  (IT 
=  B(ir 
=  E  (IT 
(C  (2,3) 
=  5  (IT 
=  3  (IT 
=  3  (IT 
=  E  (IT 


D  *C(1,1) 

,1)  *C  (2,2)  +3  (IT 
,1)  *C  (3,3)  +c  (IT 
,1)  *C  (4,4)  4E  (IT 
,1)  *(C  (1 ,2)  +C  (2 
,1)  *(C  (1,3)+C  (3 
,1)*(C(2,3)4C(3 
♦C  (2,  1)  ) 

,1)  *(T(1,4)4C  (4 
,1)  *(C  (2,4)  +  C  (4 
♦C  (2,1)  ) 

,1)  *(C  (3,4)*C  (4 
+C  (3,1)  ) 

,2)  *C  (2,2) 

, 2)  *C  (3,3)  ♦  B  ( 
, 2)  *  (C  (2 , 3)  +C  (3 
,2)  *C  (4,4)  *9  (IT 
,4)  *c  (2,2)  *3(IT 
,2)  *{C  (3,4)  +C  (4 
♦C  (3,2)  ) 

,3)  *C  (3,3) 

,3)  *C  (4,4)  4E  (IT 
,4)  *C  (3,3)  +9  (IT 
,4) *C  (4,4) 


,2)  *(C(1,2)  +C (2 , 2) ) 

,3)*  (C  (1  , 3)  +C  (3,1)) 

,4)  *  (C  (1,4)  +  C  (4,1)) 

,1))  ♦  3  (IT ,  2)  *C  (1 ,  1) 

,1))  *3(17, 3)  *C  (1,1) 

,2)  )  43(17,2)  *  (C  (1,3)  4C(3,1)  )  4 

,1)  )  43(17,4)  *C  ( 1 ,  1) 

,2)  )  +  5  (1 7,2)  *(0(1,4)  4C  (4,1))  4 

,3)  )  43(17,3)  *(C(1, 4)  +C  (4 , 1 )  )  4 


IT, 3)*  (C  (2,3)  4C  (3,2)  ) 

,2))  +B(II,?)*C(2,2) 

,4)  *  [C  (2,4)  +C  (4,2)  ) 

,2)  4  (C  (2,4)  ♦  C  (4,2)  ) 

,3)  )  +B  (17,3)  *  (C  (2 ,4)  4C  (4 ,2)  )  ♦ 


,4)*(C(3,4)4C(4,3)) 
,3)*(C(3,4)  +  C  (4,3)  ) 


CONTINUE 


W“7IZ  (6,1050)  (  (K  (IA  ,0A )  ,  I A  =  1  , 8 )  ,  0A=  1 , 28) 

FCStrZT  ('  1 ,  8  S 1  4 . 4) 


N 1  =  1 

AC  |1,M) 

= 

Cl 

AC  (2,N1) 

S 

Q2 

AC  ( 3  , N 1 ) 

s 

C  3 

AC  («, SI) 

= 

C  4 

AC  (5  ,N1) 

05 

AC  ( 6 ,1M) 

= 

Q6 

AC  (7, SI) 

s 

Q7 

AC  (6,VD 

- 

09 

CONTINUE 

DC  160  I Z 1  =  1 , 9 
V  (171)  =  AC  (IZ1  ,N1) 
CONTINUE 


NCTI:  HA  VF  TO  ADD  NON  IISEA3  ECUATIONS 


9. 


*/ 

■ij 


S’ 


ft 


NCTfc: 


ONLY  iiV^T'Y  TENTH  VALUE  IS  PLOTTED 


N3  =  n 

DC  2  1C  it 2=1, N1,  10 
113  =  N3  ♦  1 
TI  (S3)  =  AT  (N?) 

TCI  ( K3 )  =A0  (1,N2) 

TC2  <  N3  )  =  AQ  ( it ,  ”  2 ) 

IC3(53)  =  AO  (3,  N2) 

XC4(F3)  =  AQ  (4 ,  K2) 

2 1 C  CC353NLS 

C 

W?  1 1  r  (6 , 1  Cl  5) 

CALI  KFL011  (TT,TQ1 ,  N  3  , 12 ,  *  DISPLACEMENT1) 

B3II2  (o,1071)  IT  ( 1 )  ,  IMA X, DLL? 

1 C 7  1  FCSM?!  (4QX,  1  DISPLACEXFNI  (1)  AS  A  FUNCTION  OF  TAU',//,20X 
1 1 1 = •  ,  1 F7 . E , 1 0 X ,  • T MAX= 1 , 1 F  7 . 4 ,  1 CX , • DILI* • , 1 F7. 4) 

SPILB  (3,1  CIS) 

CALI  W PLOT  1  (TT,TQ2, 5?, 12, ‘DISPLACEMENT’) 

SPIT?  (6,1075)  TI  (1)  /"MAX,  DILI 

1C7S  FCPEA1  (40  X ,  '  DI  SPLACS^Sit'1’  (2 )  AS  A  FUNCTION  OF  TAU',//,20X 
1 '1=  '  ,1  F7.4 ,107,  'TMAX=‘  ,1F7.  4,1CX,  •  E£LT=  '  ,  1  ?7 . 4) 

Split  (6,1  CIS) 

CALI  WFL3T1  (IT  ,TQ3  ,  i’3 , 1  2  ,  •  DTSFLAC  EM  ENT  • ) 

KFI1I  (6 , 1 C  eG)  TT  (1)  ,TM? X,CrLT 

1 08C  F  C  3  i! A 1  (UCX, ' DISPLACEMENT  (3)  AS  A  FUNCTION  0?  TAU',//,20X 
1 •!= • ,1f7.4 ,10X, ,1F7.4,1CX,,DELI=', 1?7. 4) 

SAITS  (c.  ,1C1S) 

cm  WFLCI1  (IT  ,^Q4  ,  N  3 , 1 2  , '  DISPLACEMENT') 

SPIII  (6 , 1 C65)  T  I  ( 1 )  , TM A  X , DELI 

1085  FOPKAT  (4GX,  *  DISPLACEMENT  (U)  AS  A  FUNCTION  C?  TAU',//,20X 
1'T  =  •  ,1F7. 4,1  OX,  *TKAX='  ,  1F7.  4, 1CX,  '  DiLC  =  '  ,  1  r  7 .  4) 

SPIIF  (6  ,  1  C  1  K)  x 

GC  TC  S98 
999  CONTINUE 
SICE 
END 

SENTFY 


.  K* 


V)\®  -oNP 


APPENDIX  V 


Mode  Shape  Integrals  of  a  Simply  Supported  Beam 


'  ^  ^  * 


k4  it4  /  sin  jirC  sin  krrS  d£ 


kV 


j=k 


J7k 


kV, 

~r~6  jk 


/  Yj  Yk  d?  =  /  sin  M  sin  IcrrC  d  g 


j  j=k 


0  j^k 


=  I  6jk 


1 


/  YjYk  d? 


kir  /  sin  jirC  cos  kir?  d£ 
o 


=  a^L  (values  of  a.t  are  tabulated  in 


(A-21) 


(A-22) 


(A-23) 

Table  4) 


FORT BAN 


0001 

0002 

00C3 

0004 

0005 

0006 

00C7 


0008 
CCC9 
QC  1C 
CC11 
0012 
0013 
0014 
GC 15 
0016 


0017 
CC  18 
OC  19 

0C20 

0021 


0022 

0023 

0024 

C025 

0026 

0027 

0C28 

0029 

0030 

0031 

0032 


IV  G  LEV  El  21 


miN 


U.'-. 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


4 

♦ 

♦ 

♦ 

♦ 

♦ 

4 


DIVEKGENCE 
SIM?L Y-SUPPCPTZC 
SPECIAL  EDITION 


4 

4 

4 

4 

4 

4 

4 

444 


FEAL 

A(4#4)  ,  P  (4,4)  ,  C  (4,4) 

PEAL 

DEL  (4,4) 

PEAL 

M  A  T  M ,  :i  U  ,  N  U  ,  M ,  F  F  ,  T 

PEAL 

LAM  (4)  ,  L 

PEAL 

DO ( 6 , 2  5  C  0 ) 

FEAL 

K  (8,28)  ,  AO  (3, 2020)  ,  AT  (2020)  ,  V 

c 

SEAL 

TQ  1  (2020)  , TT  (2020)  ,TC2  (2020) 

c 

PEAD 

(5,1000)  ( A 0  ( I U ,  1 )  ,  2 U  =  1 , 8 ) 

1000 

F  G  R  M  A 

T  (3F10.4) 

PE  AB 

(5,  1005)  ((A  (I1,J1),  01  =  1,4), 

1005 

FORMA 

I  (4F10.4) 

PEAD 

(5t 1 01 0)  M,MF,T,E 

1C  1C 

FORMA 

T  (4  E 1 0 .  3) 

PEAD 

(5,1011)  VF,L,Z1 

1  Cl  1 

FORMA 

T  (F12.4,  2F10. 4) 

N  3  =  0 


c 

c 

ETTA  =  VF*S03m (MF*L**2./ (E*T) ) 

MU  =  MF/(flF  +  M) 

NO  =  7.04*Z1*S9?.T  (1  .-MU) 

C 

WRITE  (6,1^15) 

1015  FORMAT  (M‘) 

C 

C  NOTE:  DELTA  FUNCTION 

C  B  (J,K)  ,C  (K,  I)  ...  ACE  DIFFF.PE'17 

C 

CO  1  10  IE  1  =  1, a 
CO  100  TD2=1,4 

EEL  (IDl  ,  IC2)  =  9.0 
3  (ID1,ID2)  =0.0 
C (ID1,ID2)  =0.0 
ICO  CONTINUE 

DEL < I D 1 ,IE1)  =  0.5 
LAN  (ID  1 )  =  3. 14159*ID1 
E(IDl,IDl)  =  . 5*L  AM  (TD1)  **2 

C  (ID  1,  IC1)  =  0.5*LAK  (ID1)**2 
11C  CONTINUE 
C 
C 


0033 


FMU  =  SOFT  (MU) 


FORTH AN 


OC  34 
0035 
0036 
0037 


0038 

0039 

0040 

0041 

0042 

0043 

0044 

0045 


C046 
0047 
0048 
0049 
CC50 
0051 
CC  52 
0053 
0C54 

0055 

0056 

0057 

0058 

0059 

00*0 

0061 

0062 

CC63 

0064 

0065 

0066 

0067 

0068 

0069 

CC70 


IV  G  LEVEL  21  **MN  DATE  =  771  )5 

C 


C 

CO  120  14=1,8 

00(14,1)  =  0.0 
EQ  (1 4 , 2)  =0.0 
1 2  C  CONTINUE 

C 

C  NOTE:  DO-LOO?  7C  FILL  IN  O'S  5  1*S 

C 

IC  =  4 

EC  140  I A=1 , 4 
EO  130  IE=1 , 23 
K  (IA  , IB)  =0.0 
13C  CONTINUE 

IC  =  IC  M 
K  (T A  ,  IC)  =  1.0 
140  CONTINUE 
C 

C  NOTE:  EMTF.TE5  WITH  P.ECCGNIZZD  ? A T V E F. S 

C 

EO  160  15=5,8 

IT  =  15-4 
EO  150  J5=1,4 

K(7.S,J5)  =  -  (DEL  (ir,J5)  *LA?.  (IT)  **U, 0+BZTA**2. 0*3  (IT, 05 
150  CONTINUE 

EO  155  J3=5,8 

J"’=  J3-4 

K  ( 1 5  , U 3 )  =  -  (2.0*DEI?*PNU+r.  (IT,  JT)  +  DEL  (IT,J?)  «»:*J) 

155  CONTINUE 
C 

K  (15  r  9)  =  3(Tt,1)*C  (1,1) 

K  (15  #10)  =  B(I?,1)*C(2,2)+B(IT,2)*(C(1,2)+C(2,2)) 

K  (15,11)  =  B(I?,1)*C(3,3)*-?(TT,?)5*(C(1,?)*C(3,1)) 

K  (15 ,12)  =  B  (I”’  ,  1 )  *C  (4,4)  +E  (IT, 4)  *(C(1»4)  +C  (*,  1)  ) 

K  (15 , 1  3)  =  9(1”,'')*  (C  (1,2)  *C  (2,1)  )  ♦  L<  (T~  ,  ?)  *C  ('!,*’) 
K(I5,14)  =  n  (IT,1)*  (C  (1,3)  +C  (3,1)  )  +2  (In,3)  *C  (1,1) 

K  (15 , 15)  =  B(I?,1)*  (C(2, 3) +C  (3,2)) +B(Ii, 2)  *(C(V?)+=  (i 
# E  (IT  ,  3)  *  (C  ( 1 , 2)  +  C (2  ,1 )  ) 

K  (15 , 1 6)  =  3(IT,  1)  *  (C  (1#4)  (4,1)  )  +3  (IT, 4)  *C  (1  ,1) 

K  (15 , 1 7)  =  3(ir#1)*|C(2,4)*c(u,2))*5(II,2)*(C(1,4)*:(4 
*B  (IT  ,4)  *  (C  (1 , 2)  +C  (2  ,1)  ) 

K  (15 , 1 8)  =  B(IT,1)  *  (C  (3,4)  +C  (4,3))  «■?  (IT  ,3)*(C(1,4)+C(4 
#E(I?,4)*(C(1,3)  +C  (3  ,1)  ) 

K  (15  ,1 9)  =  B  (IT ,  2)  ♦C  (2,2) 

FT  (15  ,20)  =  B(IT,2)*C(3,3)  ♦  B  (I?  ,  3)  *  (C  (2, 3) +C  (3,  2)  ) 

K  (15  #21)  =  B(I7,2)  *  (C(2,3)  ♦C(J,*>)  )+1  (I  3)  *•?  (2 ,1) 

K  (15  ,22)  =  B(I?,2)*C(4,4)+?  (IT, 4)  *  ( C  (2,4)  ♦  C  (4,2)  ) 

K  (15 , 23)  =  3(1", 4)  *  c  (2,2)  ♦  F  ( IT ,  2)  *  ( C  (2,4)  +  C  (4,2)  ) 

K  (15,24)  =  ?  (I", 2)  *  (C  (3,4)  +C  (4,3)  )  +  3  (I.’  ,3)  ♦  (C  (2,a)  *Z  (4 
*8(11,4)  *(C  (2, 3)  +C  (3,2)  ) 

K  (15,25)  =  3  (IT  ,  3)  *C(3,3) 


C071 


L-v-:'  ICETRAN  IV  G  1 1 V  5  I  21 


MAIN 


n  t  ~  *  “ 

k*  A  .<  — 


s* 


e 

0C72 

K  (15  » 26)  =  B  (1^,3)  *C{4#4)  *?  (IT,  d)  * 

B 

C073 

K  (15,27)  =  B(TT#4)*C(3#  ?)♦&(!?, 3)  » 

0074 

r 

K  (15  #23)  =  P  ( I f  4)  *C(4r4) 

[£ 

0C75 

16C 

CONTINUE 

C 

v- 

C 

2  CONTINUATION  CARDS  JUST  IS  ' 

S 

0076 

165 

r* 

CONTINUE 

«_*7- 

0077 

HEAD  (5,1012,END=9?9)  T,DELT,T*!A 

» ST* 

* 

0078 

1 C  1 2 

FOFVAT  (3F10.4) 

0079 

WHITE  (6/  5554)  T,  D  ELT  ,  TMAX 

0030 

5554 

FORMAT  (  '  •  ,  3F20.  4) 

& 

CC31 

p 

Ml  =  1 

0082 

170 

CONTINUE 

r 

V 

0083 

C 

CO  160  IZ1=1,8 

CC84 

V(IZl)  »  AQ(IZ1,N1) 

ooes 

180 

CONTINUE 

c 

c 

c 

c 


NOTE:  H A  V 2  TO  A  CD  SON  i: 


:at.  equations 


•  •  t 

0086 

v  (9)  = 

?! 

Q  (1  ,N1)  **3 

0087 

V(1G) 

= 

A0(1»K1)  *AQ  (2  ,  V  1 )  *  *  2 

0088 

V  <  1 1 ) 

= 

AC  (1, Nl) ’"AO  (3 ,  Nl )  **  2 

•I*' 

.  CO  89 

V  (12 ) 

— 

AQ  ( 1  f  N 1 )  *AQ  (4,  VI)  **2 

•m 

COSO 

v  (1 3) 

S 

AC  (1,S1)  ’**2  ’♦AQ  (2,  VI) 

C0S1 

V  (14) 

- 

AC  (1  #  N 1 )  * AO  (3, SI) 

0092 

V  (1 5 ) 

= 

A  C  ( 1  /  N 1 )  *  A  Q  (2  »  N 1 )  *  A  Q  ( 3  /  :7 1 ) 

».'* 

0C93 

V  (1b) 

= 

AQ(1f  SI)  **2*f-.0  ( 4  ,  N 1 ) 

V 

0C94 

v  (17) 

= 

AC  (1  f  Nl)  *AQ  (2,ri)*Ap  (4  f  ?•  1 ) 

0095 

V  (1 8 ) 

= 

«C  (If  Si)  *A  Q  (?,>71)*AQ  (4, VI) 

m 

0096 

V  (19) 

= 

AC  ( 2 , N 1 ) **3 

0097 

V  (20) 

= 

AO  ( 2 ,  N 1 )  *  AO  (3,  VI)  **2 

* ,  • 

0098 

v  (21 ) 

= 

AO (2,M1) **2 *10  (3, vi) 

v 

0C99 

V  (22) 

= 

AC  (2, SI)  *  A.  0  <4,N1)**2 

01  CO 

V  (2  3) 

= 

AC  (2fvi)  **2*'r.Q  (4, VI) 

b’r 

0101 

V  (24) 

— 

ao  (2, si)  *ao  (.I, rl)*  aq  ( 4 , i ) 

— 

01C2 

V  (25) 

— 

AC  (3, vi) **3 

id 

0103 

V(26) 

s 

AC  (3,  Nl)  *  AC  (4  /  V 1 )  *  *  2 

•  . 

01C4 

V  (27 ) 

= 

AC  (3,  HI)  **2*?  Q  (4,  Ml) 

t . 

>v 

»v 

01C5 

C 

V  (28) 

AC  ( 4 ,  f 1 1 ) **3 

1 

C1C6 

C 

A  T  ( N  1 ) 

— 

T 

pa 

0107 

NNN  = 

Nl 

+  2 

V 

C1C8 

C 

Nl  =  N 1  ♦  1 

c 

c 

c 


710  5 

D  ) 

1)  ) 


♦ 

♦ 


MATF.IX  ■VULTIFIICAIICN 


♦ 

+ 


FOE! RAN  IV  G  LEVEL  21 


ii: 


MAIN 


D  AT  E 


771  ?r 


01C9 
0110 
cm 
0112 
0113 
0114 
Cl  15 
0116 

0117 


0118 

0119 


0120 

0121 

0122 

0123 

0124 

0125 

0126 

0127 

0128 

0129 

0130 

0131 

0132 

0133 


B 


0134 
0135 
0136 
0137 
0138 
01  39 


CO  200  IZ2= 1 , 9 
EUfcMY  =  0. 

EO  190  123=1,26 

MATM  =  K  (IZ2,I7.3)  *V  (123)  . 

LUMMY  =  EUMMY  +  MAT  F. 

19C  CONTINUE 

EO  (1 22 , NN N)  =  DUMMY 

AQ  (1 22 , N  1 )  =  DEL?*  (1.916  7*  EC  (1 22  ,  NNN)  -1.33331{D0(IZ 
#0.41 67*DQ  (IZ2 ,  NSN-2 ) )  +  AQ  (122,1:1-  1) 

20C  CONTINUE 
C 

1  =  7+DZLT 

IF  (IMAX.GT.T)  GO  TC  170 
C 

c 

EO  210  N2=1,t;1,10 

N3  =  N3+ 1 

IT  (N  3)  =  AT(N2) 

TQ1  (N3)  =  A C  ( 1  ,  N  2 ) 

102  (N3)  =  AC (2, N2) 

2 1 C  CONTINUE 
C 

EO  220  IS= 1,3 

EO  (I  S ,  2)  =  DO  (IS,  NN  K- 1 ) 

EO  (I  S,  1 )  =  CO  (TS,NMN-2) 

AO  (IS,  1)  =  AQ(IS,N1) 

220  CONTINUE* 

N3  =  N3-1 
GO  TO  165 
C 

959  CONTINUE 
C 

c 


C 

C  NOTE:  ONLY  EVERY  TENTH  VALUE  IS  PLOTTED 

C 
C 

WRITE  (6,5555)  3  3,71(1)  ,TT  (N2) 

5555  FORMAT  (•  •  ,  16,  2F1C-.  u) 

CALL  EFIPS1  PT, TCI,  S3, 1,10,0) 

CALL  DFIPSl (TT,TQ2, S3 , 1 ,  *0 .0) 

WRITE  (6,1015) 

STCP 


I] 


111 


